Method for preparing a radiation therapy plan

ABSTRACT

A method for determining a radiation treatment plan for a radiotherapy system providing multiple individual rays of intensity modulated radiation iteratively optimized the fluence of an initial set of such rays by a function that requires knowledge of only the prescribed dose and the dose resulting from the particular ray fluences. In this way, the need to store individual dose distributions of each ray are eliminated.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional application No.60/095/535, filed Aug. 6, 1998.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENTBACKGROUND OF THE INVENTION

The present invention relates generally to radiation therapy planningfor the treatment of tumors and suitable for radiation therapy machinesproviding independent intensity modulated narrow beams of radiation.

Radiation therapy involves the treatment of tumorous tissue with highenergy radiation according to a treatment plan. The treatment plancontrols the radiation's placement and dose level so that the tumoroustissue receives a sufficient dose of radiation while the radiation tosurrounding and adjacent non-tumorous tissue is minimal.

Intensity modulated radiation therapy (IMRT) treats a patient withmultiple rays of radiation each of which may be independently controlledin intensity and/or energy. The rays are directed from different anglesabout the patient and combine to provide a desired dose pattern.Typically, the radiation source consists of either high-energy X-rays,electrons from certain linear accelerators, or gamma rays from highlyfocused radioisotopes such as Co⁶⁰.

Methods of producing intensity modulated rays of radiation are wellknown in the art and include the stop and shoot method, (Xia, P.,Verhey, L. J., “Multileaf Collimation Leaf Sequencing Algorithm forIntensity Modulated Beams with Multiple Static Segments,” MedicalPhysics, 25:1424-34 (1998)), the sliding window technique (Bortfeld, etal., “Realization and Verification of Three-Dimensional ConformalRadiotherapy With Modulated Fields,” Int'l J. Radiat. Oncol. Biol.Phys., 30:899-908 (1994)), intensity modulated arc therapy, (Yu, C. X.,“Intensity-Modulated Arc Therapy With Dynamic Multileaf Collimation: AnAlternative to Tomotherapy,” Physics in Medicine & Biology, 40:1435-49(1995)), and sequential (axial) tomotherapy, (Carol, et al., “TheField-Matching Problem as it Applies to the Peacock Three DimensionalConformal System for Intensity Modulation,” Int'l J. Radiat. Oncol.Biol. Phys., 34:183-87 (1996)).

One highly accurate IMRT method uses a planar fan beam which orbits thepatient in the plane of the beam to treat a single slice of the patientat a time. Prior to reaching the patient, the fan beam is passed througha multileaf collimator (MLC) consisting of a series of opaque leaves. Asthe radiation source rotates around the patient, the tungsten leavesmove into and out of the radiation beam modulating the intensity ofindividual rays of the fan beam.

An intensity value for each ray of the fan beam at each angle of the fanbeam about the patient and for each slice of the patient is defined by atreatment sinogram. The treatment sinogram is prepared by a physicianbased on a dose map indicating the amount of radiation dose and itslocation throughout the patient.

Preparation of a treatment sinogram from a dose map is extremelycomplicated. Examples include simulated annealing (Langer M. And MorrillS., “A Comparison of Mixed Integer Programming and Fast SimulatedAnnealing For Optimized Beam Weights in Radiation Therapy,” MedicalPhysics, 23:957-64 (1996)), linear programming (Langer M. and Leong J.,“Optimization of Beam Weights Under Dose-Volume Restrictions, Int'l. J.Radiat. Oncol. Biol. Phys., 13:1225-60 (1987)), non-linear programming(Bortfeld et al., “Methods of Image Reconstruction From ProjectionsApplied to Conformal Radiotherapy” Phys. Med. Biol., 35:1423-34 (1990)),mixed-integer programming (Langer M. And Morrill S., “A Comparison ofMixed Integer Programing and Fast Simulated Annealing For Optimized BeamWeights in Radiation Therapy,” Medical Physics, 23:957-64 (1996)), anditerative filtered backprojection (Holmes et al., “An Iterative FilteredBackprojection Inverse Treatment Planning Algorithm for Tomotherapy,”Int'l. J. Radiat. Oncol. Biol. Phys., 32:1215-1225 (1995)). Anothermethod is the “Dynamically Penalized Likelihood” method suggested byLlacer and described in U.S. Pat. No. 5,602,892.

Many of these methods place severe burdens on computer memory. Forexample, in tomotherapy applications, a medium sized radiation treatmentplan will often involve storing intensities of over 91,000 rays ofradiation. Tracking the dose provided by these rays may require storageof more than 2.7×10¹¹ dose elements.

BRIEF SUMMARY OF THE INVENTION

The present invention provides a method and apparatus for generatingtreatment sinograms from dose maps.

More specifically, the present invention provides a method foroptimizing a radiation treatment plan for a radiotherapy machineproviding independently controlled radiation along a plurality of rays jdirected towards a patient to deliver dose D_(i) ^(d)=d_(ij)w_(j)tovoxels i. In a first step, a prescribed total dose D_(i) ^(p) at thevoxels i in a treatment area is received from a physician and a fluencew_(j) value is assigned to each ray j. An actual total dose D_(i) ^(d)produced at each voxel i with the assigned fluence values w_(j) is thencalculated. The fluence values w_(j) are then modified according to anupdate function of the prescribed dose D_(i) ^(p) and the actual doseD_(i) ^(d) without reference to the dose per energy fluence, d_(ij),delivered to each voxel by the given ray j. Finally the modified fluencevalues w_(j) are used to control the radiotherapy machine.

Thus it is one object of the invention to provide a method ofdetermining fluence values of multiple rays used in a radiation therapysession without the need to store partial dose values for each ray.

In one embodiment, the update function may be a ratio of the prescribeddose D_(i) ^(p) and the actual dose D_(i) ^(d) for each voxel ireceiving radiation from the given ray j or for example:$w_{j}^{({k + 1})} = {w_{j}^{k}\frac{\sum\limits_{i}\quad {aD}_{i}^{p}}{\sum\limits_{i}\quad {aD}_{i}^{d^{k}}}}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification of the fluence of the rays and a is apredetermined approximation of dose per magnitude of energy fluence,d_(ij).

Thus it is another object of the invention to provide a computationallysimple method of modifying ray fluences such as may be rapidly executedon an electronic computer. By using an approximation of dose per energyfluence, d_(ij), or dose per any magnitude related to energy fluence,the above described problems of storing and calculating partial dose areavoided.

In an alternative embodiment, the update function may be a ratio of theprescribed dose D_(i) ^(p) and the actual dose D_(i) ^(d) for each voxeli receiving radiation from the given ray j or for example:$w_{j}^{({k + 1})} = {w_{j}^{k}\frac{\left( {\prod\limits_{i}^{n}\quad D_{i}^{p}} \right)^{\frac{1}{n}}}{\left( {\prod\limits_{i}^{n}\quad D_{i}^{d^{k}}} \right)^{\frac{1}{n}}\quad}}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification of step (d).

Thus it is another object of the invention to provide a function formodifying fluences of the rays to converge to produce the desired dosehaving no partial dose d_(ij) term.

The foregoing and other objects and advantages of the invention willappear from the following description. In the description, reference ismade to the accompanying drawings which form a part hereof and in whichthere is shown by way of illustration a preferred embodiment of theinvention. Such embodiment does not necessary represent the full scopeof the invention, however, and reference must be made to the claimsherein for interpreting the scope of the invention.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

FIG. 1 is a perspective view of the shutter system assembly used in thepresent invention showing the shutter leaves and their associatedactuators;

FIG. 2 is a cross section of the shutter system of FIG. 1 along line 2—2showing the trapezoidal aspect of each shutter leaf for a radiation fanbeam of radiation, and the guide rails for supporting the shutter leaveswhen they move;

FIG. 3 is a block diagram showing the elements of a radiation therapymachine incorporating a conventional CT scanner and the shutter systemof the present invention and including a computer suitable forcontrolling that shutter system per the present invention.

FIG. 4 is a simplified representation of the gantry of the radiationtherapy machine of FIG. 3 showing variables used in the calculation of apatient model;

FIG. 5 is a flow diagram for the process of optimization of the rayfluence values per the present invention.

FIG. 6 includes two illustrations of a DVH-based system useful inguiding the optmization. FIG. 6a is an illustration of a DVH-basedpenalty system disclosed by Bortfeld et al, “Clinically RelevantIntensity Modulation Optimization Using Physical Criteria,” presented atXII International Conference on the Use of Computers in RadiationTherapy, Salt Lake City, Utah, USA, 1997 (unpublished), wherein theshaded region corresponds to the zone being penalized. FIG. 6b is ageneralization of the Bortfeld DVH penalty. Each region considered has adifferent weight in the penalization scheme.

FIG. 7a is a dose distribution of a treatment plant wherein a weight of0.95 was assigned to the treatment volume and a weight of 0.05 wasassigned to the sensitive areas.

FIG. 7b is the cumulative dose volume histogram corresponding to thedose distribution of 7 a.

FIG. 8a is a dose distribution of a treatment plan including a DVHspecification requiring a penalty to be added if more than 15% of thesensitive area exceeded a dose of 0.4.

FIG. 8b is the cumulative dose volume histogram corresponding to thedose distribution of 8 a.

FIG. 9a is a dose distribution of a treatment plan wherein a DVH basedpenalty was applied if more than 25% of the sensitive area exceeded adose of 0.1.

FIG. 9b is a the illustrated solution to the objective functioncorresponding to the dose distribution of 9 a where the solid line isthe sum of the squared differences between delivered and prescribeddoses over all pixels in the tumor and in the sensitive area. The dashline is the value of the same computation over all tumor pixels and onlythose pixels in the sensitive area that are penalized.

FIG. 10a is a dose distribution of a prostate treatment plan where thecentrally located target includes the prostate and seminal vesicles.Above the prostate is the bladder and below is the rectum. The dash lineis the 95% isodose line.

FIG. 10b is the cumulated dose volume histogram corresponding to thedose distribution of 10 a. The two specifications for the rectum areshown with square DVH and the two specifications for the bladder areindicated with diamonds on the DVH.

DETAILED DESCRIPTION OF THE INVENTION Radiotherapy Equipment

Whereas the present invention finds use with any radiation therapymachine capable of irradiating a patient at multiple angles with a largenumber of fluence controlled narrow beams of radiation in the preferredembodiment, the invention makes use of a multi-leaf collimator-typesystem.

Referring to FIG. 1, such a radiation therapy machine 10 includes aradiation source 12 producing a generally conical radiation beam 14′emanating from a focal spot 18 and directed toward a patient 17 (notshown in FIG. 1). The conical radiation beam 14′ is collimated by arectangular opaque mask 16 constructed of a set of rectangular shuttersystem blades to form a generally planar radiation fan beam 14 centeredabout a radiation fan beam plane 20.

A shutter system 22 is centered in the radiation fan beam 14 and aboutthe radiation fan beam plane 20 prior to the radiation beam beingreceived by the patient 17, and includes a plurality of adjacenttrapezoidal leaves 30 which together form an arc of constant radiusabout the focal spot 18. Each leaf is constructed of a denseradio-opaque material such as lead, tungsten, cerium, tantalum orrelated alloy.

The leaves 30 are held in sleeves 24 so that each leaf 30 may slidecompletely within its corresponding sleeve 24 to block the ray 28associated with that sleeve 24. Preferably, the leaves 30 of the shuttersystem 22 subtend the entire radiation fan beam to divide the radiationfan beam into a set of adjacent slab-like rays 28 at offset angles f.When the leaf 30 blocks its corresponding ray 28, it is referred to asbeing in the closed state. The sleeves 24 are of ample length to permiteach leaf 30 to slide out of the path of the radiation fan beam so as toleave its corresponding ray 28 completely unobstructed and yet to stillbe guided by the sleeve 24. In this nonlocking position, a leaf isreferred to as being in the “open state”.

Each leaf 30 may move rapidly between its open and closed state by meansof a primary relay-like electromagnetic actuator 32 connected to theleaf 30 by a slider member 34. The fluence passed by the ray 28 may becontrolled by changing the duty cycle of the movement of the leaf thatis the ratio of time between which it is in the open state as opposed tothe close state.

Referring to FIG. 2, the leaves 30 are supported and guided within thesleeves 24 by guide tongues 36 fitted into grooves 38 cut along theedges of the leaves 30. The grooves 38 allow the guide tongues 36 toslidably retain the leaves 30 within the sleeves 24 during motionbetween the open and closed states.

Referring now to FIG. 3, the radiation source 12 is mounted on a gantry44, the latter rotating within the radiation fan beam plane 20 about acenter of rotation 45 in the patient 17 so that the radiation fan beam14 may irradiate a slice of the patient 17 from a variety of gantryangles θ. The radiation source 12 is controlled by a radiation controlmodule 48 which turns the radiation beam 14 on or off under the controlof a computer 51.

A shutter system control 52 directed by a timer generating desiredposition signals provides electrical excitation to each electromagnet tocontrol, separately, the actuators 32 to move each of the leaves 30 inand out of its corresponding sleeve 24 and ray 38 (see also FIG. 1). Theshutter system control 52 moves the leaves 30 of the shutter system 22rapidly between their open and closed states to either fully attenuateor provide no attenuation to each ray 28. Gradations in the fluence ofeach ray, as needed for the fluence profile, are obtained by adjustingthe relative duration during which each leaf 30 is in the closedposition compared to the relative duration during which each leaf 30 isin the open position for each gantry angle.

The ratio between the closed and open states or the “duty cycle” foreach leaf 30 affects the total energy passed by a given leaf 30 at eachgantry angle and thus controls the average fluence of each ray 28. Theability to control the average fluence at each gantry angle permitsaccurate control of the dose provided by the radiation beam 14 throughthe irradiated volume of the patient 17 by therapy planning methods tobe described below. The shutter system control 52 also connects withcomputer 51 to allow program control of the shutter system 22 to bedescribed.

An optional tomographic imaging system 11 employing an x-ray source 46and an opposed detector array 50 may be advantageously mounted on thesame gantry 44 as the radiation source 12 to produce a tomographic orslice image of the irradiated slice of the patient 17 prior to radiationtherapy for planing purposes or during treatment. Alternatively, suchtomographic imaging may be performed on a separate machine and theslices aligned according to fiducial points on the patient 17.

A gantry control module 54 provides the signals necessary to rotate thegantry 44 and hence to change the position of the radiation source 12and the gantry angle q of the radiation fan beam 14 for the radiationtherapy, as well as for the computer tomography x-ray source 46 anddetector array 50 also attached to gantry 44. Gantry control module 54connects with computer 51 so that the gantry may be rotated undercomputer control and also to provide the computer 51 with a signalindicating the gantry angle q to assist in that control.

Control modules for the tomographic imaging system 11 include: x-raycontrol module 56 for turning on and off the x-ray source 46 and dataacquisition system 58 for receiving data from the detector array 50 inorder to construct a topographic image.

An image reconstructor 60 typically comprising a high speed arrayprocessor or the like receives the data from the data acquisition system58 in order to assist in “reconstructing” a tomographic treatment imagefrom such data according to methods well known in the art. The imagereconstructor 60 may also use post-patient radiation detector signals 59to produce a tomographic absorption image to be used for verificationand future therapy planning purposes as described in more detail below.

A terminal 62 comprising a keyboard and display unit 63 allows anoperator to input programs and data to the computer 51 and control theradiation therapy machine 10 and the tomographic imaging system 11 andto display images provided by the image reconstructor 60 on display unit63.

A mass storage system 64, being either a magnetic disk unit or tapedrive, allows the storage of data collected by the tomographic imagingsystem 11 and the post-patient radiation detector 53 for later use.Computer programs for operating the radiation therapy machine 10 willgenerally be stored in mass storage system 64 and laded into theinternal memory of the computer 51 for rapid processing during use ofthe radiation therapy machine 11.

The radiation source 12 may be a linear accelerator excited in pulsedmode with the pulses in synchrony with the digital to analog converterof the data acquisition system 58 so as a set of views may be obtainedduring shutter opening and closing. If each projection of radiation at agiven gantry angle q during radiotherapy is one second, the pulse rateof linear accelerator may be two hundred times per second providingreal-time motion study of movement of the leaves 30 based on thechanging fluence exiting the leaf and entering the patient 17.

During operation of the radiation therapy machine 11, the shutter systemcontrol 52 receives from the computer 51 a treatment sinogram comprisinga fluence profile for each gantry angle θ. The treatment sinogramdescribes the intensity or fluence of each ray 28 of the radiation beam14 that is desired for each gantry angle θ at a given position of thepatient support table (not shown) as translated through the radiationbeam 14.

Referring now to FIG. 4, a shutter system provides control of a totalnumber J of 10 rays 28 identified by index variable j=1 to J. Each ray28 generated by the shutter system 22 passes through the patient 17along ray center line 66 to be detected by post-patient radiationdetector 53 having detector elements.

Treatment Planning

Referring to FIG. 5, the generation of the optimal radiotherapytreatment plan according to the present invention begins with theidentification of a prescribed dose map D_(i) ^(p) providing the amountof dose desired at different voxels i within a slice as indicated byprocess block 100. Typically these different voxels i are grouped intoareas that will include one or more areas of tumorous tissue where highdose is required and one or more areas of sensitive tissue where thedose must be limited to below a predetermined value.

The prescribed dose map D_(i) ^(p) is stored within the memory of thecomputer as an array of elements, each element holding one digitalvalue. The method for entering the dose map D_(i) ^(p) may includedisplaying the tomographic view of the patient on the display of theterminal and manually tracing around the tumorous area using atrack-ball or similar input device, as is well understood in the art.Standard area-filling algorithms may be used to transfer the dose valuesassigned to each trace region to the appropriate element in the array ofmemory representing the desired dose map. Each element of the dose mapD_(i) ^(p) defines the dose desired at one voxel i within a slice of apatient.

A fluence value w_(j) of each ray j of each beam at each gantry angle θthat will produce the desired dose at each voxel i must then bedetermined as indicated by process block 102. This process is one ofiteration; an arbitrary initial fluence value w_(j) for the rays j isselected which is then modified repeatedly until optimized values areobtained.

The closer the initial fluences w_(j) selected for the rays j are to thefinal values, the faster the optimization can be completed. For thisreason, in one embodiment of the present invention, a library of priorradiotherapy treatment plans is screened to select a treatment plan fortreating a patient having a similar arrangement of tumorous tissue andsensitive tissues. The similarity between the patient, the previoustreatment plan and the current plan will provide initial fluence valuesw_(j) for the rays which are a close approximation to the rays necessaryfor the current radiotherapy application. The library may consist ofseveral different treatment plans stored within a data storage system,such as a computer, and have a catalog of various treatment volumes ofdifferent shapes and sizes.

As represented by process block 104, the delivered dose D_(i) ^(d) thatwould be provided by the initial ray fluences w_(j) is next determinedby conventional techniques. As taught in U.S. Pat. No. 5,317,616 issuedMay 31, 1994, hereby incorporated by reference, a determination ofTerma, total energy released per unit mass may be determined along eachray based on the ray's fluence and the properties of the patient. TheTerma for a given voxel may be accumulated for each ray and each beamangle and then the total Terma for each voxel convolved with aprecomputed scatter kemel(s) to determine dose at that voxel. Thekemel(s) may represent the scatter over the range of a beam angle fromdifferent beam angles and thus in one convolution operation provide thedose calculation for all beam angles. The kemel(s) may be computed byconventional techniques such as Monte Carlo simulation. The convolutionof the Terma with the scatter kemel(s) provides an accurate account oflateral scatter which is of particular importance in cases such as headand neck or tangential-field breast radiotherapy where the irradiatedvolume is small.

Generally, the Terma of each ray is not saved nor is the partial dosedelivered to a voxel by a single ray saved, thus providing substantialmemory savings.

At process block 106, the delivered dose D_(i) ^(d) calculated atprocess block 104 is compared to the prescribed dose D_(i) ^(p) enteredat process block 100 and each ray's fluence adjusted by an updatefunction relating to a ratio of a function of the prescribed dose D_(i)^(p) over a function of the actual dose D_(i) ^(d) for each voxel ireceiving radiation from the given ray j.

In a first embodiment, the update function is a ratio of the geometricmeans for the prescribed dose D_(i) ^(p) and the actual dose D_(i) ^(d)for each voxel i receiving radiation from the given ray j, and may beillustrated as follows: $\begin{matrix}{w_{j}^{({k + 1})} = \frac{{w_{j}^{k}\left( {\prod\limits_{i}^{n}\quad D_{i}^{p}} \right)}^{\frac{1}{n}}}{\left( {\prod\limits_{i}^{n}\quad D_{i}^{d^{k}}} \right)^{\frac{1}{n}}\quad}} & (1)\end{matrix}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification and n is the total number of voxels. As can beseen from inspection of equation (1), only total dose values for thevoxels are required and the partial doses contributed by the particularrays j are not needed and thus need not be stored as noted above.

It can be shown analytically that this first ratio update method whenapplied repeatedly (by repeating process blocks 104 and 106 using ineach iteration of process block 104 the modified fluence values from theprevious process block 106), that an objective function O_(j)({overscore(w)}) tends to be optimized: $\begin{matrix}{{O_{j}\left( \overset{\_}{w} \right)} = {\sum\limits_{i}\quad {\frac{1}{d_{ij}}\left\{ {{D_{i}^{d}\left\lbrack {{\ln \left( \frac{D_{i}^{p}}{D_{i}^{d}} \right)} + 1} \right\rbrack} - D_{i}^{p}} \right\}}}} & (2)\end{matrix}$

which to a first order approximation is: $\begin{matrix}{{O_{j}^{(f)}\left( \overset{\_}{w} \right)} = {\sum\limits_{i}\quad {\frac{1}{d_{ij}}\left\{ \frac{\left( {D_{i}^{p} - D_{i}^{d}} \right)^{2}}{D_{i}^{p} + D_{i}^{d}} \right\}}}} & (3)\end{matrix}$

Alternatively, in a second embodiment, the update function for modifyingthe beam weights may be a ratio of the sum at the prescribed dose D_(i)^(p) and the actual dose D_(i) ^(d) for each voxel i receiving radiationfrom the given ray j, and may be illustrated as follows: $\begin{matrix}{w_{j}^{({k + 1})} = {w_{j}^{k}{\sum\limits_{i}\quad \frac{{aD}_{i}^{p}}{\sum\limits_{i}\quad {aD}_{i}^{d^{k}}}}}} & (4)\end{matrix}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification and α is a predetermined approximation of thedose per energy fluence (d_(ij)), or dose per any magnitude related toenergy fluence, of the given ray j being modified. Alternatively α maybe a non-constant central axis depth dose stored and then looked up toserve as an approximation for d_(ij). By not storing actual values ofd_(ij), the memory requirements are still significantly reduced. In theupdate factor, the inclusion of d_(ij) would normally serve to place thegreatest importance on those voxels receiving the highest dose. Theapproximation used may influence the rate of the conversion rate of thealgorithm, but the full dose distribution determined per iteration willmaintain the accuracy of a dose computation performed using theconvolution/superposition technique.

It can be shown analytically that when this second update method isapplied repeatedly per process block 108, (by repeating process blocks104 and 106 using in each iteration of process block 104 the modifiedfluence values from the previous process block 106), that the followingobjective function O_(j)({overscore (w)}) tends to reach optimization:$\begin{matrix}{{O\left( \overset{\_}{w} \right)} = {\sum\limits_{i}\quad \left( {D_{i}^{p} - D_{i}^{d}} \right)^{n}}} & (5)\end{matrix}$

where n is an exponent having a value of 2. In a similar approach,O_(j)({overscore (w)}) may be optimized using n having value of n>2.

This equation minimizes a sum of the magnitude of the difference betweenthe delivered doses and the prescribed doses. The convex nature of thisobjective function dictates that any local minimum is also the globalminimum. With a convex objective function such as this, the use ofstochastic optimization techniques is unwarranted.

The updating method can be further modified to make the objectivefunction more robust. Specifically, the update function can be modifiedso as to apply weighting factors to each region of the patient, per thefollowing equation: $\begin{matrix}{w_{j}^{({k + 1})} = {w_{j}^{k}\left( \frac{{\sum\limits_{{i\quad \in \quad T}\quad}\quad {C_{T}d_{ij}D_{i}^{p}}} + {\sum\limits_{i\quad \in \quad R}\quad {C_{R}d_{ij}D_{i}^{p}}}}{{\sum\limits_{i\quad \in \quad T}\quad {C_{T}d_{ij}D_{i}^{d^{(k)}}}} + {\sum\limits_{i\quad \in \quad R}\quad {C_{R}d_{ij}D_{i}^{d^{(k)}}}}} \right)}} & (6)\end{matrix}$

In this equation, C_(T) is a weighting factor assigned to a tumor area,and C_(R) is a weighting factor assigned to a sensitive area. T denotesthe tumor volume and R indicates the sensitive area. As before thevalues of d_(ij) may be approximated by a constant value α or by valueslooked up in a table approximating d_(ij).

In its application, the penalty for overdosing a voxel in the tumorvolume can be set equal to the penalty for underdosing the same voxel.It is straightforward, however, to implement weighting factors thatplace a greater emphasis on either underdosage or overdosage, thusproducing a more clinically acceptable result.

The use of weighting factors is also applicable to sensitive structures.One possibility includes optimization where underdosed voxels areassigned a weight of zero. As a result, the voxels in the sensitiveareas are only penalized if they receive a dose greater than theassigned tolerance dose.

In another embodiment, the flexibility of the iterative technique isfurther improved by considering a cumulative dose volume histogram (DVH)for each treatment volume. For a particularly sensitive structure, theuser can specify a point on the DVH that indicates both the dose limit(D_(max)) and a fraction of the sensitive structure (V_(max)) that ispermitted to exceed that limit. One possible implementation of dosevolume considerations can be based upon a technique developed byBortfeld et al, “Clinically Relevant Intensity Modulation OptimizationUsing Physical Criteria,” presented at XII International Conference onthe Use of Computers in Radiation Therapy, Salt Lake City, Utah, USA,1997, (unpublished). With a DVH based penalty, one can obtain both auniform target dose and a clinically acceptable dose distribution in thesensitive areas.

The DVH based penalty guides the optimization, but its specification isnot an absolute constraint. A weighting factor may also be added to eachDVH specification thereby increasing the penalty for a violation. Byincreasing the relative weighting factor assigned to a penalty, oneeffectively raises the importance of meeting the DVH specification.

A DVH-based penalty is particularly useful with organs that are parallelin nature. This is because with parallel coordinates, the oncologist isoften willing to sacrifice a portion of the organ to obtain a favorabledose distribution in the tumor.

The present optimization technique accounts for the DVH based penaltyand the computation of the update factors. Previously all of the voxelsin the sensitive areas were assigned a tolerance dose. Dose volumeconsiderations, however, only require the inclusion of a select numberof sensitive area voxels for optimization.

According to this embodiment, a voxel in the sensitive structure ispenalized if it receives a dose between D_(max) and D′. D′ is thecurrent dose at which V_(max) is exceeded. This is illustrated in FIGS.6a and 6 b. The penalized voxels represent the voxels of the sensitiveareas receiving the smallest excess dose above D_(max), and arepenalized because they require the smallest reduction in dose in orderto satisfy the DVH specification. Accordingly, the subset of penalizedvoxels will change with each iteration.

The penalty can be added based upon any criteria. For example, it islikely that a practitioner may choose to add a penalty if more than acertain percent of the region at risk exceeds a specified dose.Likewise, the penalty could be added to the objective function if acertain condition was not met.

Under this embodiment, the algorithm determines, once per iteration, ifthe DVH specification has been fulfilled. If the specification has notbeen met, a penalty is added to the objective function. The penalty isapplied to voxels in the RAR with the smallest excess dose aboveD^(plim). Referring to FIG. 6a, the shaded region corresponds to thesevoxels. The voxels are chosen because they require the smallest changein dose so as to meet the DVH specification. In this embodiment,equation (6) may be rewritten as: $\begin{matrix}{w^{k +} = {w^{k}\left( \frac{{\sum\limits_{i\quad \in \quad T}\quad {C_{T}d_{ij}D_{i}^{p}}} + {\sum\limits_{i\quad \in \quad R}\quad {\lambda_{i}^{DVH}C_{R}d_{ij}D_{i}^{p}}}}{{\sum\limits_{i\quad \in \quad T}\quad {C_{T}d_{ij}D_{i}^{d^{(k)}}}} + {\sum\limits_{i\quad \in \quad R}\quad {\lambda_{i}^{DVH}C_{R}d_{ij}D_{i}^{d^{(k)}}}}} \right)}} & (7)\end{matrix}$

where λ_(i) ^(DVH) serves as the DVH penalty. In the above example, theDVH penalty was applied to voxels located in the shaded region of theDVH shown in FIG. 6a.

A more generalized DVH penalty is also possible. For this approach, theDVH is divided in a series of dose regions. Each region has its ownpenalty value, λ_(i) ^(DVH), used to modify the DVH according to adesired plan. A typical shape of a DVH penalty applied according to thisoptimization method is illustrated in FIG. 6b. In this case, theoptimization process is dominated by the larger λ_(i) ^(DVH) values. Thestep function shown in FIG. 6b is a representation of the pattern ofweights that can be applied, and the regions where they are applied. Theordinate, however, does not represent the actual values.

The DVH based penalty does not provide a hard constraint, but isintended to only guide the optimization. A weighting factor can be addedto each DVH specification thereby increasing the penalty for aviolation. By increasing the relative weighting factors assigned to apenalty, one effectively raises the importance of meeting the DVHspecification.

One primary advantage to the methods and apparatus of the presentinvention is that they provide the capability of performing large scaledose optimizations while minimizing the memory requirements of thechosen computer. The methods are also flexible, robust and capable ofenhancement through the addition of weighting factors assigned to eachregion of the patient, or through the addition of dose volumeconsiderations. Because of their flexibility, the present inventionfurther benefits from its ability to work efficiently in conjunctionwith the convolution/superposition based dose computation.

In the methods described above, the update factor may be calculated byupdating only the voxels located in the primary path of given ray j.This approach will ultimately result in quicker optimization planningfor complicated radiotherapy treatments such as those used intomotherapy.

EXAMPLE 1

A radiation treatment plan was optimized for an inverted U-shapedtreatment volume surrounding a rectangular sensitive area. The U-shapedtreatment volume was cut out of a 5 cm by 5 cm square, and the sensitivearea was placed in the concavity of the U.

The update factor of Equation 6 was utilized to include a weightingfactor for the treatment volume (C_(T)) and the sensitive area (C_(R))of the patient. In this equation, C_(T) and C_(R) were set at 0.95 and0.05, respectively.

FIGS. 7a and 7 b present the results from this simulation. The use ofthe weighting factors resulted in a significant improvement in thetarget dose distribution as compared with the results obtained withoutweighting factors. It was observed that the dose distribution in thetreatment volume was improved by delivering a higher dose to a largervolume of the sensitive structure. By increasing the dose to thesensitive structure, the 90% isodose line was expanded to closely matchthe border of the target.

EXAMPLE 2

A radiation treatment plan was optimized according to the methods of thepresent invention by considering a cumulative dose volume histogram(DVH). The cumulative DVH provided a DVH-based penalty which wasaccounted for in the computation of the update factor during theoptimization process. The update factors were modified to include apenalty if a specified voxel in the sensitive structure received a dosebetween D_(max) and D*. D* was defined as the current dose of whichV_(max) was exceeded.

The characterization of the penalized voxels is illustrated in FIG. 6.The penalized voxels represent the sensitive voxels of the areareceiving the smallest excess dose above D_(max). These particularvoxels were penalized because they required the smallest reduction indose in order to satisfy the DVH specification. Accordingly, the subsetof penalized voxels changes with each iteration.

FIGS. 8a and 8 b present results from an optimization process whichutilized a dose volume specification in connection with Equation 5discussed above. For the inverted U-shaped geometry, a penalty was addedif more than 15% of the region at risk exceeded a dose of 0.4. Asillustrated in FIG. 8a, the 90% isodose line closely matches theboundary of the treatment volume.

EXAMPLE 3

FIG. 9 represents the results of a treatment optimization simulationinvolving a U-shaped treatment volume and a DVH-based penalty system. Inthis simulation, a penalty was added if more than 25% of the sensitiveareas exceeded a dose of 0.1.

FIG. 9b represents the objective function value over the course of theoptimization. The solid line depicts the value of the sum of the squareddifferences between the prescribed and the actual doses over the entiretreatment volume and the sensitive areas. The dash line is the actualobjective function that is minimized when a DVH-based penalty isemployed. Specifically, it is the sum of the squared differences betweenthe delivered dose and the prescribed dose over all voxels in thetreatment volume, plus the sum of the squared differences between thedelivered dose and the dose limit of the penalized voxels. Note thatboth of these functions decreased in value with each successiveiteration.

EXAMPLE 4

DVH specifications were also tested on a simulated prostate treatmentplan. In this case, the prostate was prescribed a dose of 80 Gy. Therectum DVH specifications were: (1) add a penalty if more than 15% ofthe rectum exceeds a dose of 25 Gy and (2) add a penalty if any voxelsare above 50 Gy. The bladder DVH specifications were: (1) add a penaltyif more than the 40% of the volume exceeds a dose of 27 Gy and (2)penalize all voxels over 54 Gy.

The results of the prostate simulation are shown in FIG. 10. Note thatthe 95% isodose line closely matches the border of the target. The fourDVH specifications are plotted in FIG. 10b.

The above description has been that of a preferred embodiment of thepresent invention, it will occur to those that practice the art thatmany modifications may be made without departing from the spirit andscope of the invention. In order to apprise the public of the variousembodiments that may fall within the scope of the invention, thefollowing claims are made.

We claim:
 1. A method for optimizing a radiation treatment plan for aradiotherapy machine providing independently controlled radiation alonga plurality of rays j directed towards a patient to deliver dose D_(i)^(d)=d_(ij)w_(j) to voxels i comprising the steps: a) identifying aprescribed total dose D_(i) ^(p) at the voxels i in a treatment area; b)assigning a fluence w_(j) value for each ray j; c) calculating an actualtotal dose D_(i) ^(d) produced at each voxel i with the assigned fluencevalues w_(j) of step (b); d) for each given ray j, modifying the fluencevalue w_(j) of step (b) according to an update function of theprescribed dose D_(i) ^(p) and the actual dose D_(i) ^(d); and e) usingthe modified fluence values w_(j) to control the radiotherapy machine.2. The method of claim 1 wherein the update function is a ratio of afunction of the prescribed dose D_(i) ^(p) in the numerator and afunction of the actual dose D_(i) ^(d) in the denominator for each voxeli receiving radiation from the given ray j.
 3. The method of claim 2wherein the update function is:$w_{j}^{({k + 1})} = {w_{j}^{k}\frac{\sum\limits_{i}\quad {aD}_{i}^{p}}{\sum\limits_{i}\quad {aD}_{i}^{d^{k}}}}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification of step (d) and α is a predeterminedapproximation of dose per magnitude of energy fluence, d_(ij).
 4. Themethod of claim 3 wherein α is a constant value.
 5. The method of claim3 wherein steps (b) through (d) are repeated in multiple iterations, theassigned fluence values of step (b) taking the modified fluence valuesof preceding step (d).
 6. The method of claim 3 wherein at step (d) onlythe voxels i along the center line of the given ray j are considered. 7.The method of claim 1 wherein the update function is a ratio of thegeometric means for the prescribed dose D_(i) ^(p) and the actual doseD_(i) ^(d) for each voxel i receiving radiation from the given ray j. 8.The method of claim 7 wherein the function is:$w_{j}^{({k + 1})} = {w_{j}^{k}\frac{\left( {\prod\limits_{i}^{n}\quad D_{i}^{p}} \right)^{\frac{1}{n}}}{\left( {\prod\limits_{i}^{n}\quad D_{i}^{d^{k}}} \right)^{\frac{1}{n}}}}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification of step (d).
 9. The method of claim 7 whereinsteps (b) through (d) are repeated in multiple iterations, the assignedfluence values of step (b) taking the modified fluence values ofpreceding step (d).
 10. The method of claim 7 wherein at step (d) onlythe voxels i along the center line of the given ray j are considered.11. The method of claim 2 wherein the update function is:$w_{j}^{({k + 1})} = {w_{j}^{k}\left( \frac{{\sum\limits_{i\quad \in \quad T}\quad {C_{T}{aD}_{i}^{p}}} + {\sum\limits_{i\quad \in \quad R}\quad {C_{R}{aD}_{i}^{p}}}}{\quad {{\sum\limits_{i\quad \in \quad T}\quad {C_{T}{aD}_{i}^{d^{(k)}}}} + {\sum\limits_{i\quad \in \quad R}\quad {C_{R}{aD}_{i}^{d^{(k)}}}}}} \right)}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification of step (d), C_(T) is a weighting factor assignedto a tumor area, C_(R) is a weighting factor assigned to a sensitivearea, and α is a predetermined approximation of dose per magnitude ofenergy fluence, d_(ij).
 12. The method of claim 11 where α is a constantvalue.
 13. The method of claim 11 wherein steps (b) through (d) arerepeated in multiple iterations, the assigned fluence values of step (b)taking the modified fluence values of preceding step (d).
 14. The methodof claim 11 wherein at step (d) only the voxels i along the center lineof the given ray j are considered.
 15. The method of claim 2 wherein theupdate function is:$w^{k + 1} = {w^{k}\left( \frac{{\sum\limits_{i\quad \in \quad T}\quad {C_{T}d_{ij}D_{i}^{p}}} + {\sum\limits_{i\quad \in \quad R}\quad {\lambda_{i}^{DVH}C_{R}d_{ij}D_{i}^{p}}}}{{\sum\limits_{i\quad \in \quad T}\quad {C_{T}d_{ij}D_{i}^{d^{(k)}}}} + {\sum\limits_{i\quad \in \quad R}\quad {\lambda_{i}^{DVH}C_{R}d_{ij}D_{i}^{d^{(k)}}}}} \right)}$

where w_(j) ^((k+1)) and w_(j) ^(k) are the fluence values before andafter the modification of step (d), C_(T) is a weighting factor assignedto a tumor area, C_(R) is a weighting factor assigned to a sensitivearea, λ_(i) ^(DVH) is a penalty value assigned to each region of thepatient, and a is a predetermined approximation of dose per magnitude ofenergy fluence, d_(ij).
 16. The method of claim 15 where α is a constantvalue.
 17. The method of claim 15 wherein steps (b) through (d) arerepeated in multiple iterations, the assigned fluence values of step (b)taking the modified fluence values of preceding step (d).
 18. The methodof claim 15 wherein at step (d) only the voxels i along the center lineof the given ray j are considered.